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repeated eigenvalues and the corresponding set of sys- 
tem eigenvectors does not span the whole state-space 
of the closed-loop system. 

The paper is organized as follows. Section 2 re- 
views the problem formulation for a linear opti- 
mal control design using direct parameter optimiza- 
tion. Analytical expressions for the evaluation of the 
quadratic cost function and its gradients with respect 
to the controller design parameters are also given in 
section 2. The current approach to evaluate these 
quantities are briefly reviewed in section 3. An alter- 
nate numerical scheme for the exponential matrix and 
convolution integrals involving exponential matrices 
is presented in section 4. Approximation methods for 
the evaluation of these matrix quantities using Pade 
series are described in details in sections 5, 6, and 7. 
A design algorithm based on these numerical solution 
schemes has been incorporated into a computer-aided 
design package described in [9], A simple design prob- 
lem to motivate the need for a numerical algorithm 
that handles degeneracies in the closed-loop system 
matrix is given in section 8. Optimal solutions are 
obtained using the proposed method and the diago- 
nalization method from [9]. The numerical algorithm 
has also been applied to the synthesis of an active 
control system for the JPL large space structure de- 
veloped under the LSCL research program [16]. Re- 
sults of this application are presented in section 9. 
Conclusions are given in section 10. 

2 Problem Formulation 

In this section, we recall the problem formulation de- 
scribed in Ly [9] for the control synthesis of a robust 
low-order controller in a linear time-invariant system. 
The system P*(s) is controlled by a constant-gain 
controller C(s) as depicted in Figure 1 where y*($) is 
the controlled output vector, t/J(s) the measurement 
output vector, ^‘(s) the disturbance input vector and 
u‘(s) the control input vector. As a consistent nota- 
tion, the superscript i is used throughout this paper 
to denote the system model at the i ih plant condi- 
tion. Note that the controller C(s) is considered to 
be fixed , i.e does not vary with the design condition. 
It is modelled as a linear time-invariant system of ar- 
bitrary order whose formulation accomodates both a 
feedforward and a feedback controller structures. Ro- 
bustness requirement in the context of our problem 
formulation is defined under the conditions that the 
control system C(s) stabilizes the plant P l (s) over a 
class of design conditions (i = 1 ,N P ). 

State equations describing the system model P'(s) 
of Figure 1 are as follows. Notice that, in the prob- 


lem formulation, we assume without loss of gener- 
alities that all the system states are initially acqui- 
escent. This assumption is not restrictive since one 
can always establish impulsive inputs together 

with the appropriate influence matrix to represent 
any state initial conditions. At the i th plant condi- 
tion, the system design model is described by equa- 
tions (l)-(3) below. 

State Equations: 

/ ±*(o = FV(0 4* CV(0 4 rV(0 m 
\* i (0) = 0 v } 

where x*(t) is a n x 1 plant state vector, u*(<) an 
m x I control vector, w'(t) an m' x 1 disturbance- 
input vector, F i an n x n state matrix, G x an n x m 
control distribution matrix and T l an n x m' input- 
disturbance distribution matrix. 

Measurement Equations: 

®J(0 = (0 + Di u u { (t) 4* D\ w w\t) (2) 

where yj(£) is a p x 1 measurement vector, H\ a 
p x n state-to-measurement distribution matrix, D\ u 
a p x m control-to-measurement distribution matrix 
and D\ w a p x m' input-disturbance-to-measurement 
distribution matrix. 

Criterion Equations: 

y\{t) = Hlx'it) 4 />■„«*(<) + ( 3 ) 

where yi(t) is a p* x 1 criterion vector, H' c a p' x n 
state-to-criterion distribution matrix, D[ u a p* x m 
control-to-criterion distribution matrix and D l cw a 
p f x m f input-disturbance-to-criterion distribution 
matrix. 

For generality, the disturbances w*(0 are modeled 
as outputs of a linear time-invariant system excited 
by either impulse inputs or white noises. In this man- 
ner, one can shape the disturbance signals to have any 
deterministic responses (e.g filtered step functions, si- 
nusoidal functions, exponentially decayed or growing 
sinusoidal functions, etc...) or to model stochastic in- 
puts with any given power spectral density functions. 
At the i th plant condition, the disturbance model is 
given by equations (4)-(5) below. 

Disturbance State Equations: 

/±U0 = ^*U0 + rU'(0 m 
= o 

where x^(f) is a n! x 1 disturbance state vector, 
if(t) a m' x 1 vector of either parameterized ran- 
dom impulses (i.e t^*( t ) = E[rf l 0 ] — 0, and 

ElrforfT] ” W 0 ), or white-noise processes with 

zero mean and covariance E[if(t)if T (r)] = W 0 S(t — 
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r). The matrix W 0 is an m' x rri diagonal positive 
semi-definite matrix, F*, an ri x n' state matrix of 
the disturbance model and an n' x m / input- 
distribution matrix. 

Disturbance Output Equations: 

w i (t) = H' w xi(t) + Dii ) i (t) (5) 

where u>*(f) is a m 7 x 1 disturbance output vector, 
HI, an m* x n 7 disturbance output matrix and D' w 
an m' X m 7 direct feedthrough distribution matrix. 

State model of the controller C(s) in Figure 1 is 
that of a linear time-invariant system described by 
equations (6)- (7) below. 

Controller State Equations: 

/ i(0 = Az(t) + By\{t) ,gs 

l *(») = 0 

where z(f) is a r x 1 controller state vector, A a 
r x r state matrix of the controller and B a r x p 
measurement-input distribution matrix. 

Control Equations: 

«•’(<) = Cz(t) + D£(t) (7) 

where u*(<) is an m x 1 feedback control vector, C 
an ? 7 i x r control-output distribution matrix and D 
an m x p direct feedthrough matrix. 

For control-law synthesis, all the elements of the 
controller state matrices can be chosen as design pa- 
rameters and some of them can be left fixed at pre- 
assigned values. In addition, if needed, linear and 
nonlinear equality or inequality constraints can be 
established among the selected design parameters in 
order to ensure a particular design structure. For con- 
venience in the derivation of the performance index 
and its gradients with respect to the controller design 
parameters, we define a matrix C Q that groups all the 
controller state matrices (^4, B, C, D) in one compact 
form as follows, 


C 0 = 


D 

B 


C 

A 


J (m+r)x (p+r) 


(8) 


Thus, specification of the matrix C 0 will completely 
define the controller state model. Obviously, for the 
case of a static output-feedback design (i.e the con- 
troller order r = 0), we simply have C 0 = D. In 
Section 8, we will formulate a control design problem 
based on the minimization of a performance index us- 
ing the controller C(s) defined in equations (6)-(7). 

To examine the entire class of H 2 - norm based con- 
trol problems and to handle the problem of sensitivity 
to plant modeling uncertainties, we define the objec- 
tive function given in equations (9) and (10). This 


formulation turns out to be versatile and well-posed 
for the setting of a nonlinear constrained optimiza- 
tion problem. However, depending on the types of 
disturbance model, that is whether the disturbance 
outputs w*(t) are responses to impulse or white-noise 
inputs, different definitions for the objective function 
are needed. They are given below. 

(a) For random impulse inputs rj^t): 

\ 2 {Wp /o' e2 “ ‘ (9) 

e[ y ?{tW e {t) + } dt ) 

The expectation operator E[—] is over the ensemble 
of the random variables tjJ, in the parameterized im- 
pulse inputs r}'(t) = rf 0 6(t). Control design problems 
formulated with the above performance index J{tj) 
are often classified under the category of determinis- 
tic control. Under this category are, for example, the 
familiar control problems of command tracking con- 
trol, disturbance rejection of unwanted but known 
external input signals, implicit and explicit model- 
following designs, // 2 -control to initial conditions and 
H °°- control to sinusoidal inputs. 

(b) For random white-noise inputs r) l (t): 


1 N * 
t=i 


yfMQVdijH 


(10) 


The expectation operator E a , (— ] is over the ensem- 
ble of the random processes defined in the input vari- 
ables 7 f(t) for a closed-loop system destabilized by 
a factor a*. The destabilization effectively adds a 
value a* to the diagonal elements of the closed-loop 
system matrix. With the given performance index, 
one can address the entire class of H 2 - norm based 
control design problems. For example, we can solve 
for the linear quadratic regulator design (LQ), linear 
quadratic gaussian (LQG) design, loop transfer re- 
covery (LTR), closed-loop transfer recovery (CLTR), 
model reduction based on a minimization of H 2 - norm 
of the error. 

Note that the performance indices given in equa- 
tions (9) and (10) are evaluated to a finite-time hori- 
zon if. The use of a finite time plays a significant role 
in the implementation of a reliable design algorithm 
for the optimum steady-state solution. It should be 
recognized that the objective function is well-defined 
regardless of whether the feedback control-law is sta- 
bilizing or not. Furthermore, a class of problems as- 
sociated with command tracking of neutrally stable 
or unstable target responses (e.g step and ramp com- 
mands, sinusoidal trajectories) are only tractable un- 
der the setting of a finite-time objective function but 
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not in the confine of a steady-state objective func- 
tion where tj — ♦ oo. In practice, steady-state results, 
whenever possible, are usually achieved when the ter- 
minal time tj is equal to five or six times the slowest 
time constant in the closed-loop system responses. 

There are other unique features, besides the con- 
cept of design to a finite terminal time </, that we 
have incorporated into the design objective function 
of equations (9)-(10). First of all, this objective func- 
tion is not the usual quadratic cost function defined 
in traditional linear optimal control problems. It is 
instead a weighted average of quadratic performance 
indices evaluated over the entire set of design condi- 
tions (i = 1 ,N P ). Different weights are assigned to 
each plant condition through the scalar variable W* 
where W' > 0. Of course, \{ N p = 1, then we re- 
cover the usual quadratic cost function for a single 
nominal design condition. The time-weighted factor 
e 2a t further allows us to impose directly a stability 
requirement for the closed-loop eigenvalues. Namely, 
when a steady-state design has been achieved and the 
optimum objective function is bounded, then closed- 
loop system eigenvalues for the controllable modes 
will have real parts less than —a'. Finally, the weight- 
ing matrices Q' and R' are symmetric and positive 
semi-definite matrices. Note that our solution ap- 
proach to the minimization of the objective function 
J{tj) based on nonlinear optimization does not re- 
quire the control weighting matrix R‘ to be positive 
definite. In fact, in some design problems such as 
command tracking and model reduction, an objective 
function representing simply the tracking or model- 
matching errors does not include the control term, 
hence R * = 0. 

In this section, we provide analytical expressions 
for the evaluation of the objective function J(tj) and 
its gradients dJ/dC 0 with respect to the controller 
matrix C„. Details of the derivation can be found in 
[9], For simple technicality, the problem formulation 
assumes that there is no implicit-loop paths within 
the feedback control system. Namely, the control in- 
put u'(t) or the measurement output yj(<) should not 
have any direct link to itself. This translates into the 
conditions that one of the product DD\ U or D' SU D 
must be zero. This is not restrictive since in practice 
either actuation or sensor dynamics would be incor- 
porated into the design models and thereby result- 
ing in a system that satisfies the above assumption; 
or one can simply reformulate an equivalent problem 
with a set of measurement outputs where D' su = 0. 
Let’s assume without loss of generality that we have 
the case of DD' iU = 0, then the closed-loop system 
can be written in the form [9] shown in figure 2 or 


simply, 

= F'V*(0 + r V(f) with **(<>) = 0 ( n ) 


where 


'](n+r+n')x( 

n+r+n') 


F*+ 

G'C 

(r*+ 

G'DH\ 


G'DD',JH' W 

B( 1+ 

A+ 

B{ I-t- 

D\ u D)Hl 

BCi u C 

D i ,uD)D', ul H' w 

0 

0 

ft 


[F ](n+r+n')x m' — 


(r* 4 - G' DD\ W )D' W 
B(\ + D\ U D)D\ W D' W 

r* 

x w 


( 12 ) 


(13) 


W]pX(n+r+n') — / 

[(I + D\ U D)H\ D\ U C D\ u D)D\ w Hi] 

[D';Um‘ = [(l + Di u D)Di w Di] (15) 

[^e’]p'x(n+r+n') = 

[//’ + D[ U DH\ D’ cu C (D' CU DD\ W + D' CW )H' W ] 

(16) 


and 


[C (i ] mx (n+r+n') = [■ DH\ C DD\ W H' W \ (17) 

With these definitions, equations (2), (3) and (7) for 
yj(0.2/c(0 and u *(0 become 

= + (18) 

y*(/) = + + U 9 ) 

u i (t) = C' i x' i (t) + DD' sw D i w i 1 '(t) (20) 

For a well-posed performance index </(*/), product 
of direct feedthrough term (D x CXi DD\ w + D' CW )D' W in 
the criterion output y l c (t) and the penalty weighting 
matrix Q* must be zero. Similarly, product of the 
direct feedthrough term DD^Dl in the control out- 
put ti’(<) and the penalty weighting matrix R l must 
also be zero. Under these circumstances, the perfor- 
mance index J{t/) in equations (9) and (10) can be 
written as, 


= \Y J W;Trace{L i {t s )VW i J' T ) (21) 

1 = 1 


where 


£’(</) = /o' e (F,+a,/) ‘ 

[H' C T Q'H' C + C iT R i 0)6^' +a ' I)Tt dt 


( 22 ) 
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In the derivation of the analytical gradients of the 
performance index J(tj ) with respect to design pa- 
rameters in Co of the controller state matrices, it is 
convenient to express the closed-loop system matrices 


in terms of C 0 

explicitly, as suggested in [9], 


F ,{ -- 

= F* + (G i 0 + T 7 C 0 D\)C 0 H' 0 

(23) 

v* = 

= r 0 + (G i 0 + T 2 C 0 D\)C 0 D i 0 

(24) 


H'J = H{ + D^CoHi 

(25) 

& 

= [DH\ C DD\ W H' W ] 

(26) 


= TyCoK 


where 

[^o](n+r+n')x (n+r+n') = 

P*o](n+r+n')x (m+r) = 
[r<i](n+r+n')x m' — 


F' 0 T*Hl 
0 0 0 
0 o F' w 


G' 0 
0 I 
0 0 

' t*dL 
0 

ri. 


[^ol(p+r)X (n+r+n') — 


HI 0 D' SW H' W 
0 I 0 


[H[] p . x (n+r+n') = [H* 0 


[■®o](p+r)X m‘ — 


D' D' 

sw VJ 


D',u 0 


L^l](p+r)X (m+r) — Q 0 

[f^lp'x (m+r) = [D cu o] 
[7)]mx (m+r) = P 0] 

0 0 

I 

p2](n+r+n')x (m+r) 


p3](p+r')x (n+r+n') — 


0 I 
0 0 

0 0 0 
0 I 0 


(27) 

(28) 

(29) 

(30) 

(31) 

(32) 

(33) 

(34) 

(35) 

(36) 

(37) 


where 


X%) = 

j-tj e (F"+al)t p/ijy«p/»T e (F'' + aI) T < fa 


( 39 ) 


£*(<;) = 

r‘/ c (F"+ai) T < (H'^QH’j + C' iT RC' i )e( F " ,+0, V i dt 
0 ‘ (40) 

M%) = s?! t 9 J'*+‘*n~K8* p QH* + 

C' iT RC' i )e( F " +al)t r' i W i r HT e( F ‘ i+a V T ‘’ dtrdt 

(41) 

From equations (21)-(22) and (38)-(41), evaluation 
of the performance index J{tj) and its gradients 
dj / dC 0 would require algorithms for efficient com- 
putation of the following two types of integrals of 
matrix exponentials: X(t) = f 0 e Ar Be Cr dr and 

M{t) = J* e A ^ v ~^Be Cv De Es ds dv. In the next 

section, we review briefly the numerical algorithm de- 
veloped in [9] for the computation of X(t) and A4(t)- 


3 Current Method for Evalu- 
ating X(t) and M(t) 

Previous methods for evaluating X(t) and A4(<) in- 
volve basically the diagonalization of the matrices A , 
C and E in the exponential functions. The pro- 
cedure requires the determination of the eigenval- 
ues and eigenvectors of these matrices. It is further 
assumed for convenience that similarity transforma- 
tions can be constructed from these eigenvectors to 
diagonalize the respective matrices. Namely, there 
exist nonsingular transformations V4 and Vq an d Ve 
such that 

A = V a A a V~ 1 , C = VcAcVc 1 , E = V E A E V£ l 

(42) 

where the matrices A^, Ac and A E &re diagonal. Un- 
der these assumptions one can express for example 
the exponential function of e AI as 

e >w = e v A A A v;'t = v^VX 1 (43) 


It has been shown in [9] that derivative of the per- 
formance index J{tj) with respect to the controller 
matrix C 0 (i.e dj /dC 0 ) can be obtained explicitly 
from the following set of equations, 

dJ/dCo = 

Ya=\ W 'p{( D ?QHc + Tj RC , ')X i (if)H' 0 T + 

(G' 0 + T 2 C 0 D\ ) t [A4'(< / )/7* t + C'{t s )Y' i W i 0 D'?]+ 
T?[M i (tf)Hi T + r(< / )r*VU-D* T ](D' 1 C 0 ) T ]} 


Usage of this decomposition in the calculation of A (<) 
is shown below. 

X(t) = [ e Ar Be Cr dr 

Jo 

= V A {j\ A * T Be AcT drj VJ 1 (44) 

where B = V^ 1 BV C - Advantage of this approach 
is based on the fact that the exponential function 
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of a diagonal matrix is also diagonal. In this case, 
time integration in X(t) can be performed directly by 
explicit integration of product of scalar exponential 
functions. The resulting numerical algorithm is quite 
accurate and efficient, provided that the transforma- 
tion matrices Va and Vc are not ill-conditionned. A 
similar procedure can also be applied to the evalua- 
tion of M(t). Complete discussion can be found in 
Appendices C and D of [9]. However, breakdown of 
this algorithm will occur when the matrices A , C 
or E become degenerate or near degenerate; a situa- 
tion that becomes eminent when we address control 
of flexible structures with densely packed modes as 
demonstrated in the design examples of sections 8 
and 9. 

Clearly, in order to have a reliable design algorithm 
for optimal low-order output-feedback control synthe- 
sis [9], one must develop a robust numerical scheme 
to evaluate matrix integrals of the form shown in X (£) 
and M(t) for the case of a degenerate system. 


4 Alternative Approaches for 
Solving X(t) and M{t) 

One rather simple approach is to evaluate 

X(t) = f e Ar Be Cr dr , (45) 

Jo 

M(t) = t [ V e Aiv - s) Be Cv De Es dsdv (46) 

Jo Jo 

directly using techniques based on numerical quadra- 
ture. Efficiency of numerical integration techniques 
is poor; especially when it requires small integration 
step size for satisfactory accuracy in the case of stiff 
system matrices A } C and E. Another possibility is 
to use some types of algebraic Lyapunov equations 
for the solution of X(t) and M{t ). For example, it 
can be easily shown that the matrix X(t) can be ob- 
tained from the solution of the following Lyapunov 
equation, 

AX(t) + X(t)C=[e AT Be CT ] t Q (47) 

Solution of equation (47) exists if A,(,4) 4- Aj(C) / 0. 
This condition will not be satisfied in general for arbi- 
trary system matrices A and C . Thus, from practical 
purposes X(t) and likewise M{t) cannot be solved 
from a scheme based on Lyapunov equations. 

Another possible approach is based on the direct 
use of exponential matrix. It is well-known [12] that 
convolution integrals involving matrix exponentials, 
as represented in the matrices X(t) and M{t), can be 


derived from the matrix exponential of an augmented 
matrix. It can be shown that the matrix X(t) can 
be derived from a product of the following matrix 
exponentials, 


X(t) = e At [ I 0] exp 


{[■ 



(48) 


Thus, computation of T(<) now involves the compu- 
tation of a matrix exponential. A simple reliable algo- 
rithm for computing the matrix exponential is given 
in section 5. 

In a similar fashion, one can express the matrix 
M{t) in terms of a submatrix of a matrix exponential. 
To see this, we start from its definition 


M(t) 


n e M«-‘)Be Cv De Ea dsdv 

_ 

J* e~ Aa | J* e Av BeCv dv j De E> 

-jM/*"**'*} 

*De Es ds 


ds 

(49) 


Let’s perform a change of integration variable v — 
t — r. We have, 


M{t) = J e_y * 5 1 J eA( * r)5eC(< 


*De E ‘ ds 


^(<-■0 II € Ar Be Cr dr 


*e ct De- E{t ~ a) d[t - s)e Et 


} 


J* e Aq | J\~ Ar Be~ Cr dr 

*e Ct De~ Eq dq e Et 

* e A (l-r) Be Ct/2 e -Cr dr 

*e ct ' 2 De E(t ^dq (50) 

Notice that part of the integrand in equation (50) 
delimited by braces can be replaced by terms involv- 
ing the exponential of an augmented matrix. This 
follows simply from results developed for the matrix 
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,¥(<). With this substitution, we obtain 


M(t) = exp {[ 0 

*e c ‘/ 2 Z?e B(, - ?) dq 


Be Ct ^ ‘ 
-C 


'} 


= / 0 V° lexp {[o 

J Je c ‘/ a D^e B »d 9 




(51) 


f 

' A Be ct ' 2 0 

1 

■ 0 ■ 

[7 0 0]exp < 

0 -C e c, / 7 D 

4 

0 

l 

o 

o 

J 

7 


In this section, we have shown that the matrices X(t) 
and M(t) can be formulated in terms of the solu- 
tions of some matrix exponentials. Their evaluation 
depends therefore strongly on the accuracy and re- 
liability of numerical methods for computing matrix 
exponential. We will present one such algorithm in 
section 5. However for computational expediency, 
special consideration must also be taken to ensure 
the efficiency of the overall scheme when the upper 
limit t is large and one of the matrices A, C or D 
is unstable. Also one must economize memory re- 
quirements associated with high dimensionality of the 
augmented matrix when computing the matrix expo- 
nential. These considerations will be elaborated in 
sections 6 and 7 where we give precise algorithms for 
the computation of the matrices X(t) and A4(i) re- 
spectively. 


5 Numerical Method for the 
Matrix Exponential 

Several numerical methods are available for the com- 
putation of the matrix exponential [11]. Among all 
these, an approximation method based on Pade se- 
ries is found to be satisfactory [12]. An important 
component in any numerical routine for matrix ex- 
ponential is the scaling of the matrix argument prior 
to the series calculation. Due to the simple result 
that e At = (e Ai f 2 ) 2 , a scale factor in terms of pow- 
ers of two (i.e 2 m ) is often used. In this scheme, one 
can recover the actual value of the original matrix 
exponential by performing m squarings on the ma- 
trix exponential of the scaled matrix. The index m 
is determined based on the desired size for the scaled 
matrix. In our algorithm, scaling is applied to the 
original matrix until its oo-norm | falls below 
1/2. 


As mentioned above, the preferred series approxi- 
mation in our computation of the matrix exponential 
is the Pade series. Let’s review some of the unique 
features associated with the Pade series for the case 
of a scalar function T{z). On its most basic terms, 
it is a rational function of z of a preselected order 
that approximates the function T{z). For a given 
choice of the order of the numerator (say N) and of 
the denominator (say M), the Taylor series represen- 
tation of this Pade series must match the power series 
representation of P(z) for the first (N+M + l) terms. 
Namely, 




E 
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(52) 


In fact, the most common form of the Pade series 
is known as the diagonal sequence where the numer- 
ator and the denominator have the same order (i.e 
M = N). While it is known that the Pade series for 
the matrix exponential (i.e T{z) = e l ) converges only 
slightly faster than the Taylor series for a scalar argu- 
ment, the improvement is more significant for matrix 
argument. In the matrix case, Pade series involves 
computation of a numerator matrix Af(At) and of a 
denominator matrix 'D(At). For a diagonal Pade se- 
ries of order N, we have 


D(M) = I + (iff). (!)!$ ■* 

+ (££;& <*>*+-■ 
. (zA .... 
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and 
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The matrix exponential is simply given by 

e At =V- l {At)M’(At) 


(53) 
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Invertibility of V(At) is usually ensured by proper 
scaling of the matrix argument At. 

Another important consideration in the Pade series 
is its length N. Assuming that the matrix At has 
been scaled such that is less than 1/2, the 

parameter N can be choosen according to [12] such 


that 


r,3-2JV 


(N\)> 


(2N)\ (2N + 1)! 


< 6 


I — 


( 56 ) 


where e is a given desired tolerance for accuracy. 



With a Pade series of N terms where N is deter- 
mined from above, the approximation can be thought 
of as the exact calculation of a matrix exponential for 
a "nearby” matrix (At + E) where E is the error ma- 
trix with ||£||oo < f|M<||oo- The relative error of the 
approximation is bounded by the following inequality, 


|| e 04«+£) _ e X!| 

II^Mloo 


< c||,4<|| 00 e'll > “ll~ (57) 


and straightforward procedure to compute the matrix 
exponential of any arbitrary matrix using the Pade 
series discussed in section 5. However, it becomes a 
nontrivial task when we try to implement an efficient 
algorithm that examines carefully the issues related 
to accuracy, speed and memory requirement. The 
basic difficulties lie in the fact that the matrix ex- 
ponential is for an augmented matrix of a particular 
form, 


Thus, reducing the oo— norm of the matrix At would 
indeed improve the numerical accuracy of the ma- 
trix exponential. It has also been shown that meth- 
ods by series approximation yield better accuracy if 
the matrix argument has been preconditioned. Addi- 
tional improvement may therefore be gained by first 
preconditioning the original matrix. Another imme- 
diate benefit of lowering the oo-norm of the matrix 
being exponentiated is that the actual scaling fac- 
tor m needed would also be smaller; thereby result- 
ing in a fewer number of matrix multiplications in 
the squaring procedure. As usual, preconditioning 
a matrix tends to bring singular values of that ma- 
trix closer together (i.e. lower the condition num- 
ber), thus avoiding situation where scaling factor is 
predominantly determined by a few large singular val- 
ues, and causing significant loss of precision related 
to the set of small singular values. The most common 
method used in the precondition of a matrix is the Os- 
borne’s method [14], which minimizes the Frobenius 
norm of that matrix (and thus indirectly lowering its 
oo-norm). However, extensive tests conducted so far 
seems to indicate that preconditioning of a matrix did 
not yield significant reduction in the oo-norm and a 
smaller scaling factor to justify the added computa- 
tional efforts incurred in the Osborne’s method. The 
procedure of preconditioning a matrix is nonetheless 
recommended from the point of view of improved ac- 
curacy (see [15] and[17] ). 

In the implementation of our design algorithm for 
optimal low-order controller synthesis [9], a value of 
e = 10" 8 has been selected requiring therefore a 4- 
terms Pade series (i.e N = 4) in the evaluation of the 
matrices X*(t), C'(t) and M*(t) of equations (39)- 
(41). Additional considerations in the implementa- 
tion of the proposed method for computing X(t) and 
M{t) are given in sections 6 and 7. 
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where A = C T = F fi + a'I according to our problem 
in equations (39) and (40) for the matrices X l (tj) and 
C'(tf) respectively. Clearly, if the system matrix A 
is stable (i.e all the eigenvalues of the matrix A have 
negative real parts) then one could easily encounter 
numerical overflow when evaluating the term e~ At 
even though the matrix integrals X(t) and C(t) of in- 
terest are perfectly well-behaved. The overflow prob- 
lem occurs most likely in the final squaring process. 
To arrive at a feasible approach in the evaluation of 
X(t ), one needs to examine in details the steps taken 
in arriving at the matrix exponential of the original 
matrix starting from that of a scaled matrix (i.e in 
the squaring process). 

Let’s assume that one has scaled the input matrix 
A by AAt where At is a reasonably small time inter- 
val given by At = t/n = t/ 2 m . Thus, we need to first 
evaluate 

exp { ^ ^ ® J Af j where At = t/n = t/2 m . 


For notation convenience, we define 
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Furthermore, let W = exp(AAt) = D *. Now we can 
write our result as follows, 


X(t) = \V n [D n - l E+ D n ~ 2 EF 

+D n ~ 3 EF 2 + • • • + EF n ~ l ] y 


(60) 


6 Detailed Algorithm for 
Computing A(t) 

As seen in the previous section, the matrix .V(£) can 
be evaluated in terms of a matrix exponential as 
shown in equation (48). Conceptually, it is a simple 


or 


X(t) = W[E+WEF+W 2 EF 2 + - • - + W n - l EF n ~ 1 ]. 

(61) 

rhe above results are produced by performing m 
D ^ ^ and taking the appropriate 


iquarings of 


0 F 


abmatrix for X(t). In our application (cf. equations 


8 



(39)-(40)), the solution would therefore involve prod- 
ucts of matrices of size 2(n + r + n'). Close examina- 
tion of equation (61) leads to the following algorithm 
involving only product of matrices of size (n + r + n ) 
with the final result achievable in m steps, 


Step 1 : 

Pi = W, 

Step 2 : 

P 2 = Pi 


Q i — 

Q2 — <2i + 
PiQiRu 


R\ — F 
R 2 = R\ 


Step m : 
Pm — Pm - 1 


Qm — Qm — l“b Rm — R m -\ 

Pm-lQ m — 1 Rm- 1, 


Finally, ,T(J) = It should be noted that one 

can ” absorb” this extra factor of W (= e A ^*) into the 
matrix Q i without any change to the above algorithm 
(i.e starting the above algorithm with Q l — WE in- 
stead). This removes the need to retain the matrix 
W throughout the computation. 

Finally, one notes that the terms P> or Ri for 
(i = l,m) may underflow and become a null matrix 
for some *; in particular when the scaling factor is 
large (i.e m large). When this situation happens, one 
can simply truncate the series calculation for X(t) 
up to the i th step in the above algorithm since all 
of the significant (and nonzero) terms have already 
been accumulated into the matrix Q,. 


7 Detailed Algorithm for 
Computing Mit) 


Here the numerical algorithm is a bit involved com- 
pared to the one given for the calculation of X{t). 
This is largely due to the increased complexity of 
the argument of the matrix exponential. Following 
the procedure described in section 6, let’s perform a 
scaling upon the input matrix A by AAt such that 
computation of the matrix quantities H, J, P,U 
and W = is well-behaved. These quantities are 
defined from the following matrix exponential, 


exp 
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Due to the possible numerical underflow in the matrix 
e Ct f 2 for large t, the matrices H and J are computed 


directly from the following definitions, 

H= e Ar Be CT dre' 0 ** (63) 

Jo 


and 


J = 



e CT De Er 


dr 


(64) 


However, the computation of M 0 in equation (62) can 
still underflow due to its explicit dependence on e c ^ 2 . 
For the calculation of the matrix M{t), ideally it can 
be obtained from m squarings of equation (62). If 
carried out in this manner, potential numerical over- 
flow is eminent since, according to our equation for 
in (41), we have A T — C — E T = F 1 -h a* I. 
Hence, if the matrix C is stable, then the matrix ex- 
ponential e~ ct = V n will become unbounded. To 
bypass this difficulty, as in the calculation for X(t) } 
one needs to conduct the squaring algorithm explic- 
itly. It can be shown that the matrix M(t) can be 
computed as 


M(t) = P n ~ l M 0 + P n ~ 2 M 0 U 

+ P n ~ 3 M oU 2 + • • ■ + PM 0 U n ~ 2 
+ M 0 U n ~ l / 65 x 

+ hw 2 j + hw 3 ju + phw 3 j k } 

4- PHW A JU + P 2 HW 4 J + * 

+ P n ~ 2 HW n J 


This formulation no longer involves the matrix V , 
The above series for M can be distinguished into 
two parts — one that contains the matrix M 0 and the 
other that does not. The terras involving Mo can be 
thought of as 


' p 

M 0 ‘ 

n 

' o ‘ 

0 
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( 66 ) 


which can be performed by m squarings. The remain- 
ing terms involving H , J, W, P, and U are computa- 
tionally intensive and are of the form 


n- 2 n— 2 

EL P i HW ?+i+3 JU j , where 2 + i + j < n. 

izz 0 ;=0 

. (67) 

This equation, owing to the restriction 2 + i + j < n, 
is not easily calculated in m(= log 2 (n)) steps. A rea- 
sonably efficient procedure for computing the final 
matrix M(t) is to merge both the easily computed 
portion given in equation (66) and the more difficult 
series in equation (67) into a sequence of in steps, as 
shown in figure 3. Due to potential numerical under- 
flow, the term W % ~ 2 is not accurately obtained from 
the product W X V 2 where V — W“ l . Indeed one 
needs to recompute the term W n+2 ~ 2 ’ at each step 
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or 


of the above algorithm. This could become the ma- 
jor drawback of our scheme even though we have used 
an efficient matrix exponentiation routine to compute 
W' requiring at most 2 * logi(i) matrix multiplies. 

If in addition W n is zero (or effectively so), restric- 
tion on the indices i and j of 2 + i+j < n in equation 
(67) becomes inconsequential; hence we can express 

n = - 0 2 Ej=o P* HW 2+i+ * JU i = 

(H + PHW + P 2 HW 2 + ■ ■ -)W 2 
*(J + WJU + W 2 JU 2 H ) 

resulting in a simpler algorithm involving the follow- 
ing three terms: 


' p 

M 0 ' 

n 

' 0 ' 

0 

U 


1 


(6) ( H + PHW + P 2 HW 2 + • • • + P n ~ 2 HW n ~ 2 ) 

(c) (J + WJU + W 2 JU 2 + ■ • • + W n ~ 2 JU n ~ 2 ) 

This algorithm can again be computed in m steps as 
seen in figure 4, but now there is no costly evaluation 
of a W k term at each step. 

Further simplification of the above algorithm can 
be achieved if we make use of the fact that we have 
A - E = C 7 (cf. equation (41)) and therefore 
U = P = W T . If, for some index j < m, W 2 ' (and 
likewise U 2 ’ and P 2 ’) is zero or nearly so, then this 
calculation for A4(t) is reduced to A i(t) = HjW 2 Jj 
since M m = 0 , Hj = H m , and Jj = J m . 

In the following sections, we compare the useful- 
ness of the proposed algorithm to the early algorithm 
presented in [9] in the design of low-order optimal 
controllers for two flexible mechanical systems . 


8 A Simple Two-Mass-Spring 
Design Problem 

Control of flexible mechanical systems has been of in- 
terest in recent years [18]. This problem provides us 
a simple design case where degeneracy in the closed- 
loop eigensystem can be easily illustrated. The prob- 
lem is to control the displacement of the second mass 
by applying a force to the first mass as shown in Fig- 
ure 5. At the start, it is simple to verify that the basic 
open-loop system has a pair of degenerate eigenval- 
ues at the origin. Equations for the dynamic model 
are given below, 

m x y\ - k(y 2 - yi) + u + w ^ 

m 2 y 2 - k(yi - y 2 ) 
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(69) 

where mi = m 2 = k x = k 2 = 1. For comparison, 
we have obtained an optimal second-order controller 
design of the form, 
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0 1 

A 2 \ A 22 ’ 
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1 


C = [Cn Cia]; D = [D u ] 


using both algorithms. The control design problem 
is to minimize the following H 2 - norm of the closed- 
loop transfer function T V2W between the disturbance 
w and the displacement y 2 of the second mass through 
the controller design parameters A 2 \ , A 22 , Cn » ^12 
and Dn. We start with the following arbitrary ini- 
tial design guess of A 2 1 = —2,^22 = — Ij^ii = 
O.Cn = 0.5, and D n =0. Both algorithms con- 
verge effectively to the same optimal design gains 
of A 2 i = —0.8571, A 22 = -0.9258, C n = 0,Ci 2 = 
-0.4535 and D n = -0.2449 and with an optimum 
value \\T y9W \\l = 7.71838215122. A summary of the 
resulting closed-loop eigenvalues is given in Table 1. 
The main difference between the two algorithms is in 
the CPU time for the overall computation. Results 
are obtained for a VAX/VMS- Workstation DEC-3500 
as follows: CPU time of 19.59sec with the algorithm 
based on diagonalization and 97.36sec using the pro- 
posed method. The increase in computational load is 
expected and constitutes the basic trade-off between 
reliability and speed of the solution algorithm. The 
proposed algorithm is more reliable and with this ad- 
vantage does take a bit longer in computational time. 
With the early algorithm of [9], one cannot initiate 
the search for an optimal compensator design with 
zero gains (i.e A 2 \ = A 22 = C\\ — C\ 2 — Dn = 0) 
because, in this case, the closed-loop system would 
have two pairs of degenerate eigenvalues at the ori- 
gin; one for the rigid-body mode and the other from 
the open-loop compensator poles. To alleviate this 
problem, it was suggested that one simply starts with 
any compensator design (stabilizing or not) that pro- 
duces initially a non-degenerate closed-loop system. 
Even with these considerations, it was found that oc- 
casionally the algorithm could break down due to the 
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presence of near degeneracies in the closed-loop sys- 
tem matrix. Thus, for a reliable design method, solu- 
tion algorithm must treat degeneracies as a common 
occurrence. This situation is more evident in the op- 
timal output-feedback control design for high-order 
structural models with closely packed flexible modes. 

Future consideration would be to develop a hybrid 
algorithm taking advantage of the computational effi- 
ciency of diagonalization when the closed-loop system 
matrix is not degenerate, and turning to the current 
algorithm when degeneracies are detected. System 
degeneracy can be easily checked from the condition 
number of the eigenvector matrix. 

In the next section, we present a design problem 
where degeneracies occur frequently and therefore it 
could pose a serious difficulty for the early design 
algorithm based on diagonalization. 


An optimal low-order controller is designed to 
dampen vibration of the antenna to external exci- 
tations. To evaluate the effectivenes of the control 
system, we perform the following test. The entire 
structure is agitated using the two boom-dish actua- 
tors for the first 6.4 seconds with an applied torque 
in the form of a square wave of 0.8 second in width 
and with an amplitude of 1 N-m. The control system 
is then activated right after the excitation has been 
removed, and responses of the excited structure at 
the sensors are examined. The design objective is to 
damp out the induced vibration as fast as possible 
without excessive use of controls. Note that the nat- 
ural responses of the structure will take about a few 
minutes to decay to zero (Figure 8). 

For practical implementation, the controller design 
is choosen to be of 6 </l order and has the following 
form, 


9 The JPL Large Space Struc- 
ture Control Design 

In control problems for large flexible mechanical sys- 
tems such as space structures, causes of eigenvalue 
degeneracies are usually more subtle in nature than 
the simple case presented in section 8 for a two- 
mass-spring system. The JPL large space structure 
has been carefully designed to simulate a lightweight, 
non-rigid and lightly damped structure in a weight- 
less environment [16]. The structure itself resembles 
a large antenna with a central boom-dish apparatus 
and an extended dish consisted of hoop wires and 12 
ribs (Figure 6). There are two torque actuators (la- 
belled HAl and H A10) on the boom and dish struc- 
ture to control the two angular degrees of freedom 
in pointing maneuver, and force actuators at four rib 
root locations (labelled RAl, RAA , RA1 and RA10) 
for vibration control. From the point of view of con- 
trol design, it is a challenging problem since the plant 
has many closely spaced modes and is of reasonably 
high order. There are a total of 30 modes in the ba- 
sic structural model. The flexible modes are lightly 
damped with damping ratios ranging from 0.007 to 
0.01. The two rigid-body modes have a damping ra- 
tio of 0.12. Our design concept is to use two available 
angular displacement sensors HS 1 and HS 10 of the 
boom-dish apparatus and the two torquers HA 1 and 
HA10 collocated with these sensors for control syn- 
thesis. With this selection, 20 of the flexible modes 
associated primarily with the rib motion become un- 
controllable and unobservable. These modes are re- 
moved by modal truncation from our plant synthesis 
model. Eigenvalues of the remaining 10 modes are 
shown in Table 2. 
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(70) 


The first two states in the controller model serve as 
roll-filters, limiting the control bandwidth to less than 
50 rad/ sec. In the design optimization, we have a to- 
tal of 28 design variables: 16 in the controller A ma- 
trix and 12 in the B matrix. The objective function 
for design optimization consists of a sum of weighted 
/f 2 -norms of physical response variables observed at 
different location of the structure. It is of the form 


j (t f ) = 

Li™ § { £ Qi E a + £ R i E ° [“>(</)] | 

*/—' 00 [.=1 j = 1 J 

(71) 

Note that the expectation operator E a \- ] is for a 
system destabilized by a factor a. Table 3 lists the 
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design variables y< and their corresponding penalty 
weightings Qi. Also given in the table are the con- 
trol design weightings Rj for the actuators H A 1 and 
HA10. Responses in the above objective function are 
evaluated to random disturbances of unit white-noise 
spectra applied simultaneously at all the hub and rib 
actuators. 

The design optimization begins with the following 
arbitrary initial guess on the controller matrices A 
and B, 

' -50 0 1 0 0 0 " 

0 -50 0 0 1 0 

0 0 0 1 0 0 

A - o 0-2-100 

0 0 0 0 0 1 

0 0 0 0 -4 -4 _ 

f 01 0 1 

0 0.1 


[ ? sj 

A destabilization factor a of 0.071 was used to ensure 
that all the closed-loop eigenvalues have a real part 
less than -0.071. The optimization fails to converge 
when a destabilization factor of greater than 0.075 
was selected. This difficulty seems to be in moving 
the modes at 1.68 Hz under this controller configu- 
ration, implying that additional degrees of freedom 
must be added to the controller structure given in 
equation (70). 

While the optimization convergence itself took 13.5 
hours on a VAX/ VMS Workstation DEC-3500, the 
proposed algorithm for the calculation of the objec- 
tive function and its gradients with respect to the de- 
sign parameters is robust and leads to well-behaved 
design convergence. The final optimal values of the 
A and B matrices are shown in Figure 7. Closed-loop 
eigenvalues are given in Table 4. Primary improve- 
ment is seen in the increased damping of two modes 
at 0.65 Hz. 

Closed-loop responses of the sensor and control 
variables corresponding to this design are shown in 
Figure 8. The controlled responses decay to zero in 
about 20sec after the excitation has been removed. 
Notice that the control torques are within the desired 
limits of 1 N-m; the results are obtained through ad- 
justment of the control design weights Rj in Table 
3. This design example demonstrates the usefulness 
of a design algorithm for robust low-order controllers 
using parameter optimization, and the accompany- 
ing improvement of solution reliability using the al- 


gorithms described in sections 6 and 7 for degenerate 
systems. 

10 Conclusions 

Numerical algorithms for computing matrix exponen- 
tials and integrals of matrix exponentials have been 
developed to handle cases where the system matrix is 
degenerate. Numerical optimization combined with 
the given algorithms for the evaluation of the cost 
function and its gradients with respect to the con- 
troller design parameters has well-behaved conver- 
gence even when the closed-loop system becomes de- 
generate. These algorithms have been incorporated 
into a computer-aided-design package for synthesiz- 
ing optimal output-feedback controllers. Reliabil- 
ity of the algorithm has been demonstrated using 
typical design problems encountered in the control 
of flexible structures. Clearly this algorithm when 
combined with a previous one based on diagonaliza- 
tion would enhance significantly the overall reliabil- 
ity of the optimal design procedure for low-order con- 
trollers, thereby providing an effective automated de- 
sign environment for multivariable control synthesis. 
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Eigenvalue 

Damping 

Freq (Hz) 

-0.2290 ± 

0.3397i 

0.559 

0.0652 

-0.1553 ± 

0.8480 i 

0.180 

0.1372 

-0.0786 ± 

1.2950i 

0.061 

0.2065 


Table 1: Closed-Loop Modes (2 Mass-Spring System) 


Eigenvalue 

Damping 

Freq (Hz) 

-0.09500 

± 

0.7860i 

0.120 

0.12600 

-0.08575 

± 

0.7093i 

0.120 

0.13704 

-0.02802 

± 

4. 0024 i 

0.007 

0.63701 

-0.02929 

± 

4.1844i 

0.007 

0.66598 

-0.07405 

± 

10.583i 

0.007 

1.68434 

-0.07405 

± 

10.583i 

0.007 

1.68434 

-0.11310 

± 

IO. 6 I 61 

0.007 

2.57123 

-0.11785 

± 

16.384i 

0.007 

2.67929 

-0.21365 

± 

30.520i 

0.007 

4.85749 

-0.21365 

± 

30.5201 

0.007 

4.85749 


Table 2: Open-Loop Modes of Antenna Structure 


Variable 

Qi 

Description 

RSI 

4100 

Rib #1 root velocity 

RS 4 

3950 

Rib #4 root velocity 

RS 7 

3975 

Rib #7 root velocity 

77510 

4050 

Rib #10 root velocity 

7751 

16500 

Hub angular velocity 

77510 

15600 

Hub angular velocity 

RSI 

1100 

Rib #1 root displacement 

RS4 

1050 

Rib #4 root displacement 

RS7 

1150 

Rib #7 root displacement 

RS10 

1025 

Rib #10 root displacement 

HS1 

3900 

Hub angular displacement 

HS10 

4100 

Hub angular displacement 

Variable 

Ri 

Description 

HAl 

41 

Hub torque actuator 

HA10 

40 

Hub torque actuator 


Table 3: Design Variables for Antenna Structure 
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Table 4: Boom-Dish-Controller Closed-Loop Modes 



Figure 1: A Typical Closed-Loop System with a Feedback/Feedforward Controller 


’ *'(<) 1 f F { + G i DH i , &C (T* + G' DD\ W )H' W ] [ F(<) ‘ 

i(0 = B(l + D\ u D)Hi A + BCl u C B(l + D\ U D)D\ W H' W z{t) 

<(0 J [ 0 0 F* J L *5.(0 . 

{V +G'DD i aw )D\ l> 

+ B(l + Di u D)Di w D' w tf(t) 

ri, 


Figure 2: State Model of the Closed-Loop System 


14 













Step 0 : 
Mo 

H 

J 

p 

U 

w 

Step 1 : 

Mi = PM 0 + M 0 U 

Hi = H 

Ji = J 

p 2 

U 2 

w 2 

+HW n J 

+PHW 

+ WJU 




Step 2 : 

Mi = P 2 Mi+MiU 2 

to 

II 

35 

J2 = Jl 

p 4 

U 4 

w 4 

+H 1 W n ~ 2 Ji 

+P 2 Hi\ V 2 

+W 2 Ji u 2 




Step 3 : 

m 3 = p 4 m 2 + m 2 u 4 

h 3 = h 2 

J 3 = Jl 

p a 

U 6 

W* 

+H 2 W n ~ 6 J 2 

+P 4 H 2 W 4 

+w 4 j 2 u 4 




Step j : 

Mj = P 2 ’ 1 Mj-\ + Mj-iU 2 ’ 1 

1 

Jj = Jj-1 

pV 

U v 

w 2 ’ 

+H j - l W n+2 ~ 2 ’ Jj-i 

+P v ' Hj-iW 2 ’-' 

+W v ~ l Jj-iU 2 ’~ l 




Step m : 

M m = P n/2 M m -i + M m .iU n ' 2 

H rn — H m _ 1 

j m — Jm — 1 

pn 

u n 

w n 

+H m -iW 2 J m -i 
M(t) = Mm 

+P n > 2 H m -iW n > 2 

+[V n ' 2 J m -iU n ' 2 





Figure 3: An m-Step Calculation of M{t) 
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Step 0 : 
M 0 

H 

J 

u 

p 

w 

Step 1 : 

Mi = PMo+MoU 

Hi = H 

Ji = J 

u 2 

P 2 

w 2 


+PHW 

+WJU 




Step 2 : 

M 2 = P 2 Mi+MiU 2 

Hi = Hi 

Ji = J i 

u 4 

p 4 

w 4 


+P 2 HiW 2 

+W 2 JiU 2 




Step 3 : 

m 3 = p 4 m 2 + m 2 u 4 

H 3 = H 2 


u 8 

p8 

w 8 


+P A HiW 4 

+ W*JiU 4 




Step j : 

Mj = P 2 ’-' Mj-\ 

Hi = Hj—i 

Jj = Jj-l 

u 2 ’ 

P 2 ’ 

w 2 ’ 

+Mj-iU 2,-1 

+P 2 l ~' Hj-iW 2 ’ ' 

+W 2l ~ l Jj~iU 2 ’~ l 




Step m : 

Mm = P n/ 2 M m -l 

H m — H m _ i 

Jm J m— 1 

u n 

pn 

w n 

+M m -lU n ' 2 
M = M m + H m W 2 J m 

+P n ' 2 H m -iW r 1/2 

+ W n ! 2 J m -lU n/2 





Figure 4: A Simplified m-Step Calculation of M(t) 



Figure 5: A Two- Mass- Spring Mass System 
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Top View 



Side View 


Boom 


Ribs" 


Figure 6: Antenna Structure 
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HA 1 Control (N-m) HS 1 Response (rad) 



5.343 -1.2310 x 10" 4 
6.2118 x 10" 4 4.783 

2.701 -8.1595 x 10" 4 
2.221 9.3152 x 10~ 4 

-0.5147 5.379 

1.614 -1.252 


Figure 7: Optimized Controller Matrices for LSCL Problem 



Time (sec) Time (sec) 


Figure 8: Open-Loop (solid curve) versus Closed-Loop (dashed) Responses 
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